function out = getData(dataOption,startYear)

% The yield curve data: end of month data
loadYields     = load('GSW_yieldsMd_Mat3Md.mat');

% Remove the one-month and the three month interest rate
loadYields.yieldsMd_Mat3Md = [loadYields.yieldsMd_Mat3Md(:,1) loadYields.yieldsMd_Mat3Md(:,4:end)];

% Only include bond yields with max maturity of 120 months
select                     = loadYields.yieldsMd_Mat3Md(1,2:end) <= 120;
loadYields.yieldsMd_Mat3Md = loadYields.yieldsMd_Mat3Md(:,[1==1,select]);
% Change from discrete to continuous compounding - data is already in continous compounding
%loadYields.yieldsOut(2:end,2:end) =  log(1+loadYields.yieldsOut(2:end,2:end)/100);

% %CRSP risk-free rates for h = 1 and h = 3, in annualized terms
load('crspRet.mat');
crspRet1 = 12*log(1+crspRet(2:end,2));
crspRet3 = log(1+crspRet(2:end,3));
crspRet3 = 4*filter(ones(3,1),1,crspRet3);


% Loading all GSW yields (only needed for excess return regressions)
tmp = load('GSW_yieldsMd_Mat1Md.mat');
% TEST
%tmpMagnus = load('GSW_yieldsMd_Mat1Md_2016.mat');
%tmpMagnus.yieldsMd_Mat1Md(2:end,2:end) = log(1+tmpMagnus.yieldsMd_Mat1Md(2:end,2:end)/100);
%max(abs(tmpMagnus.yieldsMd_Mat1Md(2:end,120)-tmp.yieldsMd_Mat1Md(2:end-12,120)/100))


% Load macro variables
macro = load('macroVariables.mat');

% Only include data from "startYear" and ending at 2016:12 (macro data does not go further)
% For yields
select1 = (str2num(datestr(loadYields.yieldsMd_Mat3Md(2:end,1),10)) >= startYear);
select2 = (str2num(datestr(loadYields.yieldsMd_Mat3Md(2:end,1),10)) <= 2016);
selectYield = select1.*select2 ==1;
loadYields.yieldsMd_Mat3Md = loadYields.yieldsMd_Mat3Md([1==1;selectYield],:);
loadYields.yieldsMd_Mat3Md(2:end,2:end) = loadYields.yieldsMd_Mat3Md(2:end,2:end)/100; %yields in base units

% For all yields
tmp.yieldsMd_Mat1Md = tmp.yieldsMd_Mat1Md([1==1;selectYield],:);
tmp.yieldsMd_Mat1Md(2:end,2:end) = tmp.yieldsMd_Mat1Md(2:end,2:end)/100; %yields in base units

% For macro (and CRSP yields)
select1 = (str2num(datestr(macro.macroDates(1:end,1),10)) >= startYear);
select2 = (str2num(datestr(macro.macroDates(1:end,1),10)) <= 2016);
selectMacro = select1.*select2 ==1;


%% Generating output
out.yieldsOut        = loadYields.yieldsMd_Mat3Md;
out.yields_h1        = tmp.yieldsMd_Mat1Md;
out.crspRiskFree_h1  = crspRet1(selectMacro,1);
out.crspRiskFree_h3  = crspRet3(selectMacro,1);
out.NBERrec          = macro.NBER(selectMacro,1);   
if dataOption == 1
    % we use PMI: threshold value is 44.5
    out.econ      = (macro.PMI(selectMacro,1)-44.5)/100;
elseif dataOption == 2
    % Loading the Cieslak and Pavolo inflation factor
    out.econ  = macro.CPI(selectMacro,1);
    out.PMI   = (macro.PMI(selectMacro,1)-44.5)/100;
    % Reduce sample to account for missings in CPI at the start
    idx            = ~isnan(out.econ);
    out.econ  = out.econ(idx,1);
    out.PMI   = out.PMI(idx,1);
    out.NBERrec = out.NBERrec(idx,1);
    out.yieldsOut = out.yieldsOut([1==1;idx],:);
    out.yields_h1 = out.yields_h1([1==1;idx],:);
elseif dataOption == 3
    % Loading output gap 
    out.econ = macro.GAP(selectMacro,1);
    out.PMI  = (macro.PMI(selectMacro,1)-44.5)/100;
elseif dataOption == 4
    % Loading the Ludvigson and Ng factor 
    out.econ = macro.LN(selectMacro,1);
    out.PMI  = (macro.PMI(selectMacro,1)-44.5)/100;
end 

end

